!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
 
     double precision function ran1(idum)

!     "Minimal" random number generator of Park and Miller with
!     Bays-Durham shuffle and added safeguards. Returns a uniform
!     random deviate between 0.0 and 1.0 (exclusive of the end points).
!     Call with idum a negative integer to initialize; thereafter do
!     not alter idum between successive deviates in a sequence. RNMX
!     should approximate the largest floating point value that is 
!     less than 1.
      ! use para1
      implicit none
        
      integer::idum     !=-2008
      
      integer ia, im, iq, ir, ntab, ndiv
      real am, eps, rnmx
      
      parameter (ia=16807, im=2147483647, am=1./im)
      parameter (iq=127773, ir=2836, ntab=32, ndiv=1+(im-1)/ntab)
      parameter (eps=1.2e-7, rnmx=1.-eps)

      integer j, k, iv(ntab), iy
      save iv, iy

      data iv /ntab*0/, iy /0/

      if (idum.le.0 .or. iy.eq.0) then
         idum = max(-idum,1)
         do j = ntab+8, 1, -1
            k = idum/iq
            idum = ia*(idum-k*iq)-ir*k
            if (idum .lt. 0) idum = idum + im
            if (j .le. ntab) iv(j) = idum
         enddo
         iy = iv(1)
      endif

      k = idum/iq
      idum = ia*(idum-k*iq) - ir*k
      if (idum.lt.0) idum = idum + im
      j = 1 + iy/ndiv
      iy = iv(j)
      iv(j) = idum
      ran1 = min(am*iy,rnmx)
      return
      end


FUNCTION cauchyul()
use para,only :deltaul,muul
implicit none
!INTEGER::idum
double precision::cauchyul
double precision,parameter::pi=3.141592653589793d0
double precision,external::ran1
integer,save::seedul=-1234

double precision::scale=1.0
 cauchyul=deltaul*dtan(pi*(ran1(seedul)-0.5))+muul
 cauchyul=cauchyul*scale

end function

FUNCTION cauchyub()
use para,only :deltaub,muub
implicit none
!INTEGER::idum
double precision::cauchyub
double precision,parameter::pi=3.141592653589793d0
double precision,external::ran1
integer,save::seedub=-4321

double precision::scale=1.0
 cauchyub=deltaub*dtan(pi*(ran1(seedub)-0.5))+muub
 cauchyub=cauchyub*scale

end function


FUNCTION cauchys()
use para,only :deltas,mus
implicit none
!INTEGER::idum
double precision::cauchys
double precision,parameter::pi=3.141592653589793d0
double precision,external::ran1
integer,save::seeds=-5678

double precision::scale=1.0
 cauchys=deltas*dtan(pi*(ran1(seeds)-0.5))+mus
 cauchys=cauchys*scale

end function

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
